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Abstract 

In this paper, we propose a moment method to numerically solve the Vlasov equa- 
tions using the framework of the NRxx method developed in [6, 8, 7] for the Boltzmann 
equation. Due to the same convection term of the Boltzmann equation and the Vlasov 
equation, it is very convenient to use the moment expansion in the NFLcx method to 
approximate the distribution function in the Vlasov equations. The moment closure 
recently presented in [5] is applied to achieve the globally hyperbolicity so that the 
local wcll-posedness of the moment system is attained. This makes our simulations 
using high order moment expansion accessible in the case of the distribution far away 
from the equilibrium which appears very often in the solution of the Vlasov equations. 
With the moment expansion of the distribution function, the acceleration in the ve- 
locity space results in an ordinary differential system of the macroscopic velocity, thus 
is easy to be handled. The numerical method we developed can keep both the mass 
and the momentum conserved. We carry out the simulations of both the Vlasov- 
Poisson equations and the Vlasov-Poisson-BGK equations to study the linear Landau 
damping. The numerical convergence is exhibited in terms of the moment number 
and the spatial grid size, respectively. The variation of discretized energy as well as 
the dependence of the recurrence time on moment order is investigated. The linear 
Landau damping is well captured for different wave numbers and collision frequencies. 
We find that the Landau damping rate linearly and monotonically converges in the 
spatial grid size. The results are in perfect agreement with the theoretic data in the 
collisionlcss case. 
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1 Introduction 

The Vlasov equation is a differential equation describing the time evolution of the 
distribution function of plasma consisting of charged particles with the long-range (for 
example, Coulomb) interaction. The equation was first suggested for the description of 
plasma by A. Vlasov in 1938 [16]. Due to the presence of long-range Coulomb interaction 
in the plasma, Vlasov started from the Vlasov equation, which is a collisionless Boltzmann 
equation, and used a self-consistent collective field created by the charged plasma particles 
to get the Vlasov-Maxwell equations [11]. The Vlasov-Poisson (V-P) equations are an 
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approximation of the Vlasov-Maxwell equations in the nonrelativistic zero-magnetic field 
limit. The V-P equations are used to describe various phenomena in plasma, in particular 
Landau damping and the distributions in a double layer plasma. The Vlasov-Pission- 
BGK (V-B) equations are also studied with a collision term presented as the BGK term. 
They are the simplest kinetic equations which correctly describe the essential features 
of collective and dissipative (entropy-producing) particle interactions in semiconductor 
plasmas [1, 2, 17, 18]. 

Due to the complex phenomena in the plasma, numerical simulation plays an important 
role in the study of the Vlasov and related equations. There are several kinds of methods 
to solve the Vlasov equations. The finite element method was proposed in [19, 20], which 
can be used to handle complicated boundary problems but inconvenient to solve the high 
dimension equations [10]. Meanwhile, the particle-in-cell (PIC) method [3] which used a 
finite number of macro-particles to approximate the plasma is easy to be implemented, and 
the method to discretize the Vlasov equations on a mesh of phase space was introduced 
to remedy the problem in PIC that the inherent artificial discrete particle noise made 
the description inaccuracy. The semi-Lagrangian method [15] and the cubic interpolated 
propagation method [12] were also used to solve the Vlasov equations. However, the first 
method destroyed the local characters due to the reconstruction and the second one was 
quite expensive for the storage of the distribution function and its derivatives. In [10], 
Filbet put forward a new method to deal with the force term, which made the scheme 
keep the mass and energy conserved. Recently, an approach based on the moment method 
has been proposed in [13], and therein the distribution function was expanded using the 
Hermite polynomials with a prescribed macroscopic velocity chosen as the expansion center 
and a prescribed thermal velocity as the scaling factor at different locations. 

In the past years, a regularized moment method was developed in [6] to numerically 
solve the Boltzmann equation. This method adopts the Hermite polynomial expansion to 
approximate the distribution function, with the basis function shifted by the local macro- 
scopic velocity and scaled by the square root of the local temperature. The approximated 
distribution function is used to directly solve the Boltzmann equation without the de- 
duction of the moment system up to arbitrary order of moments. The method therein 
was further explored as the NRxx method [8, 7] by introducing the regularization term 
using asymptotic expansion in term of the mean free path. Recently, a new regularization 
method [5, 4] was derived with the guarantee that the regularized moment system is glob- 
ally hyperbolic. Due to the locally well-posedness provided by the global hyperbolicity, 
it is eventually accessible that approximating the distribution function far away from the 
equilibrium distribution by the stable simulation using large number of moments. Inspired 
by this progress, we in this paper develop the NRxx method to study the Vlasov equations, 
which is similar to the Boltzmann equation in the convection term, while the distribution 
function is much farther away from the equilibrium state than the gas flows, due to the 
long-range Coulomb interactions. 

Here we are focusing on the V-P and V-B equations. With the moment expansion 
in the Vlasov equations, the convection part is smoothly handled by the original NRxx 
method for the Boltzmann equation. We discretize the regularized term given in [4] di- 
rectly using the central difference scheme. Our deduction shows that the electric potential 
brings us an acceleration on the macroscopic velocity, thus it is very convenient to be 
numerically integrated. Currently, the collision term under our consideration is the simple 
BGK model to avoid distraction. The approximated electric potential is obtained by a 
three-point central scheme, and the numerical method keeps both the mass and the mo 
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mentum conserved. By the numerical resolution study, it is exhibited that our method 
is numerically converged by the comparison of Landau damping rates obtained using dif- 
ferent spatial grid size and number of moments. With the increasing of the number of 
moments, the recurrence time is almost linearly related to the square root of the moment 
expansion order. The discretized energy of our method only varies slightly in comparison 
with the overall energy. Since the Landau damping rate is affected by the wave number 
and the collision frequency, different wave numbers and collision frequencies are exten- 
sively studied. The wave numbers ranging from 0.2 to 0.5 are simulated, and the collision 
frequencies are taken as 0.0 (the collisionless case), 0.01 and 0.05. The results show that 
our method can capture the linear Landau damping very well and the damping rates ob- 
tained is in quantitative agreement with the theoretic data. We find with surprise that 
the numerical Landau damping rate is linearly and monotonically converged in term of 
the spatial grid size. Particularly, the numerical damping rate converges to the theoretic 
data perfectly in the collisionless case. This observation inspires us to predict the Landau 
damping rates by the extrapolation of our numerical damping rates with presence of the 
collision term. 

The layout of this paper is as follows: in Section 2, the regularized moment system is 
deduced for the Vlasov-BGK equations. In Section 3, we present the detailed procedure 
of the numerical method. In Section 4, the numerical examples including the numerical 
resolution study and the linear Landau damping simulation with different parameters are 
presented. Some concluding remarks are given in the last section. 



2 Regularized Moment System 

Let f(t,x,v), which depends on time t, position x € O C M 3 and the microscopic 
velocity v G M 3 , be the distribution function depicting the motion of particles. It is 
governed by the V-B equations 

+ v ■ V x / + F(t, x, v) ■ V v f = v(J M ~ /), (2.1) 

where v(t,x) denotes the collision frequency, and F(t,x,v) is the electric force produced 
by the self-consistent electric filed E(t,x): 

F(t, x, v) = ±E(t, x), E(t, x) = -V x <f>(t, x), -A x (f) = q^, (2.2) 
m eo 

where 4>(t, x) is the electric potential produced by the particles; p, q, m and eo stand for the 
density, the single charge, the mass of one particle and the electric constant respectively; 
Jm is the Maxwellian defined as 

p(t,x) ( \v -u(t,x)\ 2 \ 

fM = (2nu th (t,x))^ 6XP {- 2u th (t,x) ) ' (2 - 3) 

where the parameter uth(t,x) is the thermal velocity [13], u is the macroscopic velocity 
and p(t, x) is the same as that in (2.2). fM is related to / by 

v )dv= [ f M (v) ( V 

M 2 / JR3 I , , 2 
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In the case of v = 0, we get the V-P equations. The relations between the macroscopic 
variables and the distribution function are deduced as 



p(t,x) = / f(t,x,v)dv, (2.5) 

p(t,x)u(t,x) = / vf(t,x,v)dv, (2.6) 
Jr 3 

p(t,x)\u(t,x)\ 2 + 3p(t,x)u th {t,x) = [ \v\ 2 f(t,x,v)dv. (2.7) 

JR 3 

The conservation of mass, momentum and total energy are all valid for the V-B and 
V-P equations. Multiplying the equation (2.1) by 1 and v, direct integration with v and 
x gives us 



d_ 
dt 
d_ 
di 



/ f(t, x,v)dxdv = 0, t£R + , (2.8) 

JR 3 xR 3 

/ vf(t,x,v)dxdv = 0, t G R + . (2.9) 

JR 3 xR 3 



Multiplying the equation (2.1) by \v\ 2 and integrating by parts, we get the conservation 
of the total energy for the system (2.1) and (2.2): 



d 
dt 



/ f{t,x,v)\v\ 2 dxdv + / \E(t,x)\ 2 dx j = 0, t G M + . (2.10) 



2.1 Hermite expansion of the distribution function 

Following the method in [6, 9], we expand the distribution function into Hermite series 



as 



/(«) = E ( 2 - n ) 

aGN 3 

where a = (a^c^a^) is a three-dimensional multi-index, and 



t= V -^. (2.12) 



The basis functions H Uth , a are defined as 

3 1 / c2- 



^ t „a(o = n ^o +1 # e <*(&)«p (-|) > ( 2 - 13 ) 

where -He^ is the Hermite polynomial 

tfe^) = (-lfexp (^) -^exp (-^) . (2.14) 

For convenience, He n (x) is taken as zero if n < 0, thus "H Uth ,a{£) is zer ° when any com- 
ponent of a is negative. With such an expansion, the Maxwellian fu can be written as 

f M (v) = pH Uth , ($), (2.15) 
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and the BGK collision term is written as 

K/m - /) = -v fofiWAZ)- ( 2 - 16 ) 

|<*|>1 

The definition of £ shows that each basis function is an exponentially decreasing func- 
tion multiplied by a multi-dimensional Hermite polynomial shifted by the local macro- 
scopic velocity u and scaled by the square root of the local thermal velocity u^. For any 
vector u' and positive number u' th , if the distribution function / is expanded as 



f(v) = £ f' a H u > h , a (t% € = "-T^, (2-17) 



then the following relations hold 

P = fd, (2.18a) 
pu = pu' + (4)J =lj2i3 , (2.18b) 

3 

p\u\ 2 + 3pu th = 2pu • u' - p\u'\ 2 + J2«hfo + 2/LJ- (2.18c) 

d=l 

In the case of u' = u and u' th = uth in (2.18) , the following relations between the coefficients 
f a can be verified: 

3 

fo = p(t,x), f ei = 0, ]T/ 2ed = 0, i = 1,2,3. (2.19) 

d=l 

Moreover, direct calculations give us the relations between the coefficients f a in (2.11) as 

3 



q i = 2h ei + Y,he d+ e i , (2.20) 

d=l 

1 3 

Pij ~ ^Pdd = (1 + %)/ ei +e r (2.21) 



where i,j = 1, 2, 3, qi and are related to / by 

<li = l [ \v-u\ 2 (vi- Ui )fdv, (2.22) 
= / (vi - Ui)(vj - Uj)fdv. (2.23) 



2.2 Moment expansion of Vlasov equation 

To get the moment system, we substitute the expansion (2.11) into (2.1) and then 
match the coefficients for the same basis functions. Taking the temporal and spatial 
derivatives directly on the basis functions H Uth ,a, the term with the expansion (2.11) 
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is expanded as 
dfc 

at 



E 



dfg Y^du d 1 du th 



=1 d=l 



E 



5 — - + «jo 1- («i + !)— 5 — 



+ 2 i u thfa-e d - ej + Ujf a - ed + (oj + l)/ Q _ ed+£j ) 



+ - 



1 duth 

2 dxi 



^ { u thfa-2e d -ej + Ujfa-2e d + («j + l)/a-2e d +e,-) 



cf=l 



the acceleration term jP • V v / is expanded as 

3 

— ^1 Fdfa-e d Hu th ,a 

and the collision term u($m — f) is expanded as 



^ XT faH-u th ,a 
\a\>l 



'Uth 

(2.24) 



(2.25) 



(2.26) 



Substituting these expansions back into (2.1), and collecting coefficients for the same basis 
functions, we get the following general moment equations with a slight rearrangement: 



df, 



/ dud 



du d 



1 / duth 



duth 



Of 



Of 



j'=i 



2 



3=1 



dxi 



-2e d 



+ E 

3 / 

+ Yl ( U th 

3=1 V 



|^ (u th f a - ed - ej + («j + l)/ a _ ed+ej ) + (uthfa^ea-ej + (otj + l)/c*-2e d +e,-) 



9/a 



5a; i 



+ u J^7 + ( a i + 1 )- 



2 



i/(l - 5(a)) f a , 



(2.27) 



5(a) 



(2.28) 



where 5(a) is defined by 

' 0, if |a| ^ 2, 
1, otherwise. 

Following the method in [7], we deduce the mass conservation in the case of a = as 

(2.29) 



21° + ( u ®A + f ^i] = 

9x j V'^j °^3'/ 



If we set a = e<j, with d = 1,2, 3, (2.27) reduces to 



/o 



c?tt 



3=1 



0, (2.30) 



3=1 
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which is simplified as 

Then we consider the case of \a\ > 2. Multiplying |t? — u\ 2 on both sides of (2.1), and 
integrating with respect to d on I 3 , we have 



, / du th J-y du th \ 2 A / dqj ^ du d \ 



j=l \ J d=l 



Remark 1. Since 



/ \v - u\ 2 ^dv = -2 f (vi - Ui)fdv = 0, z = 1, 2, 3, 

Jr3 dvi J R3 



(2.33) 



the acceleration term does not appear in (2.32). 

Substituting (2.31) and (2.32) into (2.27), we eliminate the temporal derivatives of v 
and uth- Then the quasi-linear form of the moment system reads: 

3 

+ £ 

L L "' 



— (u th f a - ed - ej + (atj + l)/ Q _ ed+ej .) + (u t / l / Q -2e d - ej + («j + lj/a-aej+ej 



3,d=l 

u(l-5(a))f a , V|a|>2 



+ E( v ^^-+^^-J+E^ + 1 ) 



(2.34) 



We collect the equations (2.29), (2.31), (2.32) and (2.34) together to obtain a moment 
system with infinite number of equations. 

2.3 Closure of the moment system 

With a truncation of (2.11), (2.34) will result in a finite moment system. Precisely, 
we let M ^ 3 be a positive integer and only the coefficients in the set M. = {f a }\ a \^M 
are considered. Let Fu(u,Uth) denotes the linear space spanned by all H Uth ,a(€ys with 
\a\ ^ M, and the expansion (2.11) is truncated as 

f( V ) « y, f* n ^ ( !L i=) • ( 2 - 35 ) 



|a|<Af 

with /(v) 6 FM(u,Uth) and f a £ M. The moment equations which contain df a /dt with 
|a| > M are disregarded in (2.34). Then, (2.29), (2.31) and (2.34) with 2 sC \a\ < M lead 
to a system with finite number of equations. 

Due to the presence of the terms with df a+ej /dxj , \a\ = M, the moment system we 
obtained is not closed yet. We rewrite (2.34) into the form below: 

^- + A a + B a = u(l-5(a))f a , (2.36) 
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where in the case of 2 < \a\ < M, 

3 3 



/0 d=l 4=1 ° X 3 4=1 V ° X 3 



+ U 



J 

3 



dfa 



2> + 1 >-a^- 

3=1 3 



Ba = 0, 

and in the case of |a| = M, 



JU d=l j=l J j=l ■> 

B a = ^2 Baj = (ay + 1)- 



<9x 



Clearly, the moments /«+ e - in term ^ Q with \a\ = M are not in the set of moments 
M. = {faja^Mi and have to be substituted by some expressions consisting of lower order 
moments to make the moment system closed. If we simply let B a = 0, the Grad-type 
system associated with the moment set M = {f a }\a\^M is obtained. It is known that 
the Grad-type system is not locally well-posed due to the lack of the global hyperbolicity, 
which results in numerical blow-up when the distribution function is far away from the 
equilibrium state. The regularization method in [4, 5] is adopted here to achieve a globally 
hyperbolic moment closure. The (M + l)-st order terms are substituted as bellow 

d fa+e 3 ( A du d 1 /A du th \\ . 

~> " (Z>-^^7 + 2 {^/^^ ) ) ' M = M, j = 1,2,3. 

(2.37) 

Let B a denote the regularization term based on the characteristic speed correction in [4, 5] 
for Id = M 



B n 



~ B »J ~ ~ ( X] fa-ea+e, ^ + £ ^ /a-2e d + ej "7^ ] • (2-38) 

d=l \d=l J d=l 3 J 



Substituting the (M + l)-st order term B a with the regularization term B a , the moment 
equations are revised as 

^ + Aa + B a = u(l-6(a))f a , (2.39) 

with B a = B a = for 2 ^ \a\ < M. If the distribution function / only depends on x\ in 
the spatial direction, we have that 

Baj = 0, for j = 2, 3. (2.40) 

Remark 2. The regularized moment system (2.39) is not able to be written into a conser- 
vation law, for the presence of the regularization term. If we let B a = with \a\ = M, 
the system changes into the conservative Grad-type system. 



S 



3 Numerical Method 

The numerical scheme for the regularized moment system (2.39) is a natural extension 
of the method in [6]. By a standard fraction step method, we split the V-B equations into 
the following parts: 



the convection step: 
the acceleration step: 



^ + v-V x f = 0, (3.1) 



df 

-± + F(t,x,v)-V v f = 0, (3.2) 

F(t,x,v) = — E(t,x), Eit.x) = -Vsc <t>(t,x), -A x (j) = q—. (3.3) 
m e 

• the collision step: 

% = K/M - /)• (3.4) 

We observe that only (2.31) contains the electric force F in the governing equations. Thus 
the governing equations of the acceleration part turn into 

d t u = F, (3.5) 

F(t, x, v) = ^E(t, x), E(t, x) = -V x <f>(t, x), -A x cj) = q^. (3.6) 
m e 

Here we restrict our study in ID spatial space. The numerical scheme adopted in the 
x-direction is the standard finite volume discretization. Suppose to be a uniform mesh 
in M, and each cell is identified by an index j. For a fixed xq G K and Ax > 0, 

T h = {Tj = x + (jAx, (j + I) Ax) :j£Z}. (3.7) 

The numerical solution which is the approximation of the distribution function / at t = t n 
is denoted as 

ft(x,v) = f?(v), x £ Tj, (3.8) 

where fj'iv) is the approximation over the cell Tj at the n-th time step and has the 
following Hermite expansion form as 

\a\^M 

3.1 The convection step 

The equation (3.1) is discretized as 

= /» + Klj(v) + K2j(v), (3.9) 

where K™ ■ is the contribution of the term A a in (2.39) without considering the acceler- 
ation, and K^j is the contribution of the term B a in (2.39). Noticing that the term A a 
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results in the conservative part in the Grad-type moment system, its contribution K™j is 
discretized in the conservative formation as 

At n r i 
Kl j {v) = -— Ff^(v) - F?M , (3.10) 

where F n 1 is the numerical flux between cell Tj and T,-+i at t n . We use the HLL scheme 

3+2 



in our numerical experiments following [8], which reads: 



^7 + i(«) 



A L , t < < A?, i , 



(3.11) 

where A L , i and A R , i are the fastest signal speeds [51 as 

3+2 J+2 



\ L +i = min{n^ - C M +i^u? hJ , u\ j+1 - C M+1 ^Ju? hJ+l }, 
Xf + i = maxj^j. + C M +iJu%~,ul j+1 + Cm+i^j+i}, 



(3.12) 



where Cm+i is the greatest zero of HeM+i(x), u\ is the first component of the macroscopic 
velocity u, and uth is the thermal velocity. The formula of the signal speed is also used 
to determine the time step At n by the CFL condition. Two points remain unclear in the 
numerical flux. The first one is how to calculate v\f™{v). This is managed according to 
the recursion relation of Hermite polynomials: 

vif?(v) = iKhj) 112 ^) + «;)] E fZoVUZ?) 

\a\^M 

_ 3.13 

\a\<M 

where £™ = (v — Uj)/ ^fu^~, £^ is the first entry of £™, and u r j,u r t l h j are the mean 

macroscopic velocity and thermal velocity in the j-th cell. Since \a + e±\ = M + 1 when 
\a\ = M, v\fj{v) no longer exists in the space Fm (it" , u™ h ■) . By simply dropping the 
terms with |a + e\\ = M + 1, we project v\f™{v) back into Fm{u™ ,u™ h since when 
|a| > M, H a {^) is orthogonal to FM{u,u t h) with respect to the inner product 



if, 9) 



^ 3 /(£MOexp(^W (3.14) 



The other point is how to add up two distribution functions lying in Fm (u™, j) and 
l> u tkj+l) respectively. The proposition in [6] is referred to solve it. 

Proposition 1. Suppose f 6 Fm{u\, u^h i) can 6e represented by 

f( v )= E /l.a%HM.a(£l)> ^1 = ~ "O/V^' ( 3 - 15 ) 

|a|<Af 
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For some u 2 G M 3 and Uth,2 > 0, {F a (T")}\aKM satisfies 

—j^T = ^2 ^ [ u th,\RF a -2e d + Wdy/U t h,lF a -e d ] , VV G [0, 1], ^ ^ 

I F«(0)=/l,a, 



where S and R are given below. And w = (u\ — «2)/\/%,2i u t h = -\/uth,i/ u th. 



2- 



R(r)= Ut \,\ v S(t) = 1-tR(t) = - L— -. (3.17) 

Kft - l)r + 1 (u^ - l)r + 1 

Lei 

= Yl F M)Ku th , 2 ,a(t2), £2 = (v - U2)/V«tM- (3-18) 

|a|<M 

T/ien <?(■«) G Fm(u 2 , 11^,2) an d d( v ) satisfies 

[ p(v)f(v)dv = [ p(v)g(v)dv, Vp(«) G P A /(M 3 ). (3.19) 

The proposition provides an algorithm to project a function in Fm(u±, uth 1) to the 
space F M (u 2 ,Uth,2)- Assuming /i G F M (ui, u th>1 ) and / 2 G F M (u 2 , u thi2 ), we can find 
an /i G FM{u2,u t h,2) as a representation of f\ G i*Af (t*i, 11^1,1) in the sense of (3.19). 
Thus, to add up fj G F M (u™,u? h j ) and G F M (w" + i, ^j+i); we first find an f j+1 G 
Fm^ 1 ] ,Uth j) as an approximation of in the sense of (3.19), and then add up the 
coefficients of fj and fj+i for the same basis function. 

The regularization term appears only in the moment equations containing df a /dt with 
\a\ = M. Therefore, only when \a\ = M, we have to calculate KJfj- We simply use the 
central difference scheme to approximate the spatial derivatives in (2.38): 



8x ~ * ^ 2 Ax 

d^tfe _ vyh^n A U th,j + 1 ~ U th,j-1 

Then we get the numerical approximation for K%j with \a\ = M, as 

Jf& = -At Y, («1 + 1) E (te^VMj + / "~ 2 2 6d+ei Vgtl&j) ?^, a 



(3.20) 




(3.21) 



3.2 The acceleration and collision step 

The acceleration step is performed by solving 

^ - *i = 0. (3.22) 



For the time step At, (3.22) is approximated as 

n+l TI-+ 



vJlf = u n 1 f>* + AtF?f, (3.23) 
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n +l l„/. n +l i „/,"+l 

j (3.25) 



where u™^ 1 '* is the first entry of the macroscopic velocity u in the j-th cell after the 
convection step at t = t n . And F^ 1 is the first entry of the electric force F in the j-th 
cell after the convection step at t = t n . In the ID spatial space, the Poisson equation (2.2) 
reduces into a second order ODE as 

F 1 (t,x,v) = ±E 1 (t,x), E 1 (t,x) = -^^-, -0 M V = -- (3.24) 
m ox €q 

The three-point central difference scheme is used to discretize the potential equation 

Ax 2 e 

where Pj +1 is the density in the j-th cell after the convection step, noticing that the density 
is not updated in the acceleration step and the collision step. The central difference is 
used to approximate E^j and F^ 1 

7 /; n + 1 _ ■>/i n + 1 

E ^ ~ 2Ax ' F U ~ m E ^ ■ (3 ' 26) 

For the BGK collision model, the moment expansion results in a simple form (2.16), 
and (3.4) changes into 

^ = -vf a , \a\ > 2. (3.27) 

It is directly integrated as 

= f^^M-^t n ), VaEN 3 , K|a|<M. (3.28) 

Remark 3. The collision step only revises the coefficients with 2 ^ \a\ ^ M, and does not 
change the macroscopic velocity u and the density p, which are decided by the coefficients 
with \a\ < 2. 

3.3 Outline of the algorithm 

The overall numerical scheme is summarized as below: 

1. Let n = and set the initial value of fj a ; 

2. Calculate At n according to the CFL condition; 

3. Integrate the convection term using (3.9); 

4- Update the macroscopic velocity using (3.26) and (3.23); 

5. Apply the collision step (3.28); 

6. Let n <— n + 1, and return to step 2. 

Remark 4. For the V-P equation, v = 0, and the collision step is simply skipped over. 
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The scheme keeps both the discretized mass and momentum conserved. Precisely, in 
the case of the periodic boundary condition, let us denote the discretized mass as 

TV 

V h {t n ) = Y j /\x f?(v)dv (3.29) 

and discretized momentum as 

TV 

M h (t n ) = Y / Ax vffdv, (3.30) 
we have the following conclusion. 

Theorem 1. The numerical solution f£(x,v) given by the algorithm above in the case of 
the periodic boundary condition satisfies that 

V h {t n ) = V h (to), M h {t n ) = M h (t ) (3.31) 

for all n > 0. 

Proof. Noticing that the mass on each cell is not modified when we apply the regularization 
term and in the acceleration step, the conservation of the mass is straight forward based 
on results in [8]. 

The momentum conservation is equivalent to verify 

M h (t n+1 ) = M h {t n ). (3.32) 

It is clear that the collision step does not change the macroscopic velocity and the density. 
Therefore we verify below that the momentum is conserved in the acceleration step and 
the convection step, respectively. 

1. We first verify that the acceleration step keeps momentum conservation. From (2.6) 
and (3.23), we have 

N 



M h {t n+1 ) = Y j Ax I vr +1 dv 

TV 

= Ax^p] 

TV 

= Ax P? +1 K +1 '* + Ff l At^ 



(3.33) 



3=1 

TV TV 

+1 



= Ax ]T P ] +1 u] +l >* + AxAt 11 F] +1 P ] 

3=1 ' 3=1 

According to (3.25) and (3.26), we have 

TV TV 



sw=i;-7 



Ax 2 V mJ 2 Ax 

TV 



3=1 3=1 

TV 

x 3 2 (^3+1^3+1 ~ ^3-1^3-1 ~ 2 ^j4>3+l ^ 3 - 34 ) 
3=1 



V>JV + ^Af+l -1pQ-ll>l+ IpQlpl - 1pNlpN+l- 
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Since we restrict the problem to the ID spatial space with the periodic boundary 
condition, we deduce by ipN+i = ipi, ipo = ipN, F2 = and F3 = that 

JV 

3=1 

and we obtain 

N 

M h {t n+l ) = AxJ2p] +1 u] +1 '*- (3-35) 
3=1 

Since p™ +1 is the density after the n-th convection step, Ax X/^=i pl +1 u^ +1 '* 1S ^ ne 
total momentum after the n-th convection step. 

2. Here we verify that the convection step does not change the total momentum. 
Thanks to (2.6) and (3.9), we have 

N N . 

Ax P n 3 +lu T 1 '* = E A * / v/; +1 '*dv 

N 

= VAx / vffdv + J^Ax / vK^dv + J^ Ax / vK^dv 

N N N 

= Ax [ vf?dv -At n J2 v(F? + i/2 ~ Ff_ l/2 )dv + £ Ax [ vK^dv 
=M h (t n ) - At n [ v(F^ +1/2 - F? /2 )dv + V Ax [ vK^dv. 

(3.36) 

Due to the periodic boundary condition, we have that FJ^ +1 , 2 = F™, 2 . Meanwhile, 
the regularization part only updates the M-th order terms which have no effect on 
the macroscopic velocity u and the density p, thus the regularization term will not 
break the momentum conservation. Precisely speaking, the basis functions of K2J 

are i0l f(v — it") j ^ u th ~j \ > w hh \a\ = M, M ^ 3, which are orthogonal to v, 

then 

vK$ d dv = 0, (3.37) 



and we obtain 

N 

Ax £>?+\ n+1 '* = M h {t n ). (3.38) 
3=1 

With (3.35) and (3.38), we conclude the total momentum conservation consequently 

M h (t n+1 ) =M h (t n ). (3.39) 
This ends the proof. □ 
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4 Numerical Examples 



We study the linear Landau damping modelled by the V-P and V-B equations with 
the periodic boundary conditions. The CFL number is always set as 0.45. The specific 
examples are from [10]. The form of the V-B equations with the periodic boundary 
conditions coupled with the normalized Poisson equation is 

df 

-±+v-V x f + E-V v f = u(f M -f), (4.1) 
E(t,x) = -V x ^(t,x), (4.2) 

-Ai/>(t,x)= f(t,x,v)dv-l. (4.3) 

Here we adopt the same initial data as in [10] 

f(0,x,v) = f M = -^=e-\ v \ 2 / 2 (l + A C os(kx)), (x, v) G (0, L) x ]R 3 , (4.4) 
V 2n 



where A is the amptitude of the perturbation, k denotes the wave number, and the periodic 
length is taken as L = 2-rr/k. The initial data and (2.35) give us 

P (0,x) = f (0,x) = (l + Acos(kx)), fa\ t=0 = 0, \a\>0. (4.5) 

What of one's interests is the evolution of the square root of the electric energy, which 
is defined as 

S h (t) = ^AxE](t). (4.6) 

i=i 

According to Landau's theory, the time evolution of the square root of £h(t) is expected 
to be exponentially decaying almost with a fixed rate jl, which is affected by the wave 
number k and the collision frequency v. For this purpose, we always plot the square root 
of £h(t) in logarithm scale in the figures in this section. The total energy is the sum of 
the electric energy and the kinetic and internal energy of the particles as 

£total(t)=£h(t)+£p(t), (4.7) 

where the kinetic and internal energy of the particles 

N 

£ p (t) = Ax^2 (/*(*)«?(*) + Pi(t)uth,i{t)) . (4.8) 

i=l 



4.1 Numerical resolution study 

We first examine the numerical convergence on different spatial grid size. The V-P 
equations in ID spatial space and ID velocity space are numerically solved. The number 
of moments is set as 80. In Figure 1, the evolution of the square root of £h(t) with the 
wave numbers k = 0.3, 0.5 on different grid size is presented. The number of spatial grids 
we used are 500, 1000, 2000, 4000 and 8000. It is clear that the square root of £ h (t) is 
damping exponentially on all different grid size. With the increasing of the grid number, 
the damping rate is decreasing monotonically. For the spatial grid size Ax, we adopt the 
least square fitting, which uses the peak value points of £h(t), to obtain the numerical 
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(a) k = 0.3 (b) k = 0.5 



Figure 1: Exponentially damping in time of the square root of £h on different spatial grids 
with the wave number k = 0.3 and 0.5. The curves in blue are the square root of Eh in 
time using logarithm scale on different spatial grid size. The slopes of the blue lines are 
the numerical damping rate 7^ by the least square fitting of the peak value points of £ h . 
The slope of the red line is the damping rate given by the theoretic data in Table 1. 




(a) k = 0.3 (b) k = 0.5 



Figure 2: The linear dependence of the numerical damping rates in spatial grid size with 
the wave number k = 0.3 and 0.5. The x-axis is the spatial grid size Ax and the y-axis 
is numerical damping rate 7^ (Ax). The line is obtained by the least square fitting of 
7^ (Ax) with Ax ranging from L/100 to L/8000. The intercept of the line on y-axis is the 
parameter 7?' and the slope of the line is the parameter 7?' 1 . 

damping rate 7^. One finds obviously in Figure 1 that both £h(t) and the damping rate 
7^; are converging while the spatial grid size Ax is going to zero. Furthermore, if we take 
the numerical damping rate 7^ as a function of the spatial grid size Ax and apply a least 
square fitting to retrieve the parameters 7^'° and 7^' 1 (see Figure 2) in the ansatz as 

ht a \ h,0 1 h.l a 
7l(Ax) = 7 L ' +7 L ' Ax, 

it is found that the fitting provides us both parameters extremely close to constants for 
all the values of k we tested ranging from 0.2 to 0.5. This indicates us that the following 
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relation 

7 £(A*)- 7 £' «A* (4.9) 

is approximately valid. The obtained parameter 7^' is regarded as the limit of the numer- 
ical damping rate when Ax is going to zero. To our surprise, the limit 7^° we obtained is 
in perfect agreement with the theoretic data in Table 1. 



Wave number k 


Theoretic data [14] 


h,0 

1l 


0.2 


-5.5 x 10~ 5 


-9.28 x 10~ 5 


0.3 


-0.0126 


-0.01260 


0.4 


-0.0661 


-0.06614 


0.5 


-0.1533 


-0.15334 



Table 1: The comparison of the limit numerical damping rates and the theoretic data. 




(a) k = 0.3 (b) k = 0.5 



Figure 3: Exponentially damping in time of the square root of £f L on different number of 
moments M with the wave number k = 0.3 and 0.5. The curves are the square root of £^ 
in time using logarithm scale using different number of moments. 

Let us turn to the study of the numerical convergence in term of the number of mo- 
ments. Again the V-P equations in ID spatial space and ID velocity space are numerically 
solved. The results using different number of moments ranging from 10 to 50 are collected 
in Figure 3 with fixed number of spatial grids. In this figure, we observe that the behavior 
of the time evolution of £h (t) for different number of moments are almost the same at the 
beginning, indicating that the numerical damping rates are very close to each other. The 
damping of the square root of £h(t) is persisting until the recurrence appears. After the 
appearance of the recurrence, the evolving of £h(t) deviates evidently from exponentially 
decaying. 

The exponential convergence rate in the number of moments is expected since the 
Hermite spectral expansion is used in the velocity space. We notice that the exponential 
convergence in term of number of moments is quite different to be observed, since the nu- 
merical error is in linear convergence due to the dominance of the spatial discretization. As 
the best try, a very fine spatial grid with 8000 points, which is about the maximal capacity 
of our current hardware, is used in the computation to suppress the spatial discretization 
error. With the fixed spatial grid size, we then take the numerical damping rate 7^ as a 
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(a) k = 0.3 (b) k = 0.5 



Figure 4: The dependence of the recurrence time on the number of moments used. The 
x-axis is the square root of the number of moments M, and the y-axis is time. The 
recurrence time for a given M is between the time of the sequential peak values of £h 
before and after it deviates from the exponentially decay (see Figure 3). In the figures, 
the vertical line segments are given by two time of the sequential peak values. 




40 45 50 55 60 65 70 75 80 ' 40 45 50 55 60 65 70 75 SO 



(a) k = 0.3 (b) k = 0.5 

Figure 5: The approximately exponential convergence of the numerical damping rates 
in the number of moments with the wave number k = 0.3 and 0.5. The x-axis is the 
number of moments M, with Mj = 40, 50, 60, 70, 80, 90 and AM = 10. The y-axis is 
log(7^(Mj) — 7^(Mj_i)), where 7^(Mj) is obtained by the least square fitting of the peak 
values of £h- The spatial grid size is fixed as Ax = L/8000. 

function of the number of moments M. In the case of an exponential convergence, the 
dependence of 7^(M) on M is as 

7£(M)= 7 r°+7r' 1 A- M , ( 4 - iQ ) 

where 7™'°, 7™' 1 and A are parameters independent of M. Let Mj be a given arithmetic 
sequence with a const AM = Mi — Mj—i. Based on the ansatz (4.10), we have that 

7 £(m_ 1 ) = 7 -'°+7r 1 A- M! - 1 - 
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By substr acting the two equations and then taking a logarithm on both sides, we have 



log ( 7 £(M) - ^(Mi-i)) = -Mi log(A) + log 7 f 1 (l - A AM ). 

It is to find that log (^y^(Mi) — 7^(Mj_i)) is approximately linear in Mj if (4.10) is valid. 
We give the plot of log (7^ (Mi) - 7 ^(Mi_i)) in M { for M { = 40, 50, 60, 70, 80, and 90 in 
Figure 5. It is clear in this figure that log (7^(Mj) — 7^(Mj_i)) is almost linearly related 
with Mi, indicating that (4.10) is valid and the numerical convergence rate in term of the 
number of moments is approximately exponential. 

It has been pointed out in [14] that the recurrence is an essential phenomenon of a 
class of numerical methods for the Vlasov equations. Such a phenomenon is mainly the 
effect of the free streaming part of the Vlasov equation. In the ID case, it is 

It is not difficult to find that for any velocity Vi, we have f(jL/vi,x,Vi) = f(0,x,Vi) for 
an arbitrary integer j when the periodic boundary condition is imposed. Therefore, if a 
numerical scheme approximates the distribution function by taking its values on vq, ■ ■ ■ , vm 
in the velocity space, then, for a time T such that most of T/(L/vi), i = 0, • • • , M are 
close to some positive integers, the discrete distribution function f(T, x, v) will be close to 
the initial setting f(0,x,v). Thus the recurrence occurs. 

Our method is not able to escape away from this trap, either, since as pointed out 
in [5], the hyperbolic moment equation is similar as a discrete velocity model with a 
shifted and scaled stencil. Based on our observation in Figure 4, it is conjectured that the 
recurrence time is proportional to the square root of the moment expansion order M. For 
a finite difference discretisation in the velocity space, the recurrence time is proportional 
to the grid size in v [14]. Noticing that the minimal distance between the zeros of M-th 
degree Hermite polynomial is around \/M, the linear dependence of the recurrence time 
on \[M consists with the analysis in [14] if we regard the moment method as a collocation 
spectral method in the velocity space with the collocation points being the zeros of Hermite 
polynomial. 

It has been proved in Theorem 1 that the mass and the momentum are conserved 
by our scheme, while the total energy is not conserved. To examine the behavior of the 
total energy of our method, we present in Figure 6 the variation of the total energy in 
time in serveral different setups. It is clear that the total energy of the numerical solution 
produced by our method is changed very slightly in the whole computation. 



4.2 Parameter study of linear Landau damping 

The Landau damping with different wave numbers and collision frequencies are studied 
in this section. It is found that the numerical phenomena here are exactly the same as the 
collisionless case presented in Section 4.1, for the wave numbers ranging from 0.2 to 0.5 
and the collision frequencies ranging from 0.01 to 0.05. We have checked the numerical 
convergence both in spatial grid size and order of moment expansion, and satisfied results 
are obtained. In all cases, the exponential damping of the square root of £h(t) is observed. 

Let us study the behavior of the solution with different wave numbers firstly. We use 
3000 spatial grids and 80 moments to compute the evolution of £h(t) for the wave numbers 
k = 0.2, 0.3, 0.4, and 0.5 with the collision frequency v = 0, and the results are in Figure 
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(c) k = 0.4 (d) k = 0.5 



Figure 6: The variation of the total energy Stotai in time in the collisionless case with the 
wave number k = 0.2, 0.3, 0.4 and 0.5. The x-axis is time and the y-axis is Stotalit) ~ 
^totality- Noticing that £totai(0) is of O(l), it is clear that the variation is very small, 
though the total energy is not conserved. 

7. It is clear the £h(t) is damping exponentially in time and the damping rate is increasing 
with greater wave number. This consists with the empirical formula given in [10]. 

In Figure 8, the evolution of the square root of £h with the wave number k = 0.2, 
0.3, 0.4 and 0.5 with the collision frequency v = 0.01 is presented. We find that £h{t) 
is also damping exponentially but with greater rates than that in the collisionless case. 
In Figure 9, the numerical results for the wave number k = 0.2, 0.3, 0.4 and 0.5 with a 
greater collision frequency v = 0.05 are plot. We find that £h(t) is damping exponentially 
with an even greater rate than that in the case of v = 0.01. This consists again with the 
empirical formula given in [10], though a different collision term was used therein. 

We numerically solve the V-B equations in ID spatial space and ID velocity space with 
the number of moments fixed as 80. The number of the spatial girds is ranging form 100 to 
8000. The least square fitting of the peak value points is again adopted to get the numerical 
damping rate of the square root of £^ for different Ax. In Figure 10, the damping rates for 
the wave numbers k = 0.3, 0.5, and the collision frequencies v = 0.01, 0.05 are plotted. It 
is clear that the numerical damping rates are all in a linearly and monotonically converging 
pattern with the spatial grid size Ax going to zero. Remember that in the collisionless 
case, the numerical damping rate is monotonically linearly converging to the theoretic 
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(a) v = 0,k = 0.2 



(b) v = 0, k = 0.3 




(c) v = 0,k = 0.4 



(d) v = 0,k = 0.5 



Figure 7: The dependence of the damping rate on the wave number k in the collisionless 
case. The curves are the evolving of the square root of £h in time using logarithm scale. 
The slope of the line is obtained by the least square fitting. 



data with the spatial grid size going to zero. Based on such similarity, it is reasonable 
to conjecture that the damping rate with the collision term as the function of the spatial 
grid size is monotonically linearly converging to the exact damping rate while the spatial 
grid size is going to zero, too. This inspires us to predict the exact damping rate for 
different collision frequencies by the extrapolation of the numerical damping rates, which 
is the same as the approach in the collisionless case. We adopt the formula (4.9) again to 
retrieve the parameters by least square fitting. The obtained parameter 7^'° is regarded 
as the prediction of the exact damping rate. In Table 2, we present 7^'° for the wave 
numbers k = 0.2, 0.3, 0.4, and 0.5 with v = 0.01 and 0.05. 



5 Concluding remarks 

The NRxx method has been extended to solve the V-P and V-B equations. The 
method is able to capture the linear Landau damping effectively. We attempt to predict 
the exact Landau damping rates with the collision term based on our observation in the 
collisionless case. 
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5 10 15 20 25 5 10 15 20 25 

(a) v = 0.01, k = 0.2 (b) i/ = 0.01, fc = 0.3 




(c) i/ = 0.01, jfc = 0.4 (d) v = 0.01, k = 0.5 



Figure 8: The dependence of the damping rate on the wave number k with collision 
frequency v = 0.01. The curves are the evolving of the square root of £^ in time using 
logarithm scale. The slope of the line is obtained by the least square fitting. 



Collision frequency v 


Wave number k 


h,0 

7l 


0.01 


0.2 


-8.3796 x 10~ 5 


0.3 


-0.01361 


0.4 


-0.06701 


0.5 


-0.15331 


0.05 


0.2 


-0.002094 


0.3 


-0.01755 


0.4 


-0.06966 


0.5 


-0.15230 



Table 2: The prediction of the linear Landau damping rates. 
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(a) v = 0.05, k = 0.2 



(b) v = 0.05, k = 0.3 
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(c) 1/ = 0.05, k = 0.4 



(d) i/ = 0.05, fc = 0.5 



Figure 9: The dependence of the damping rate on the wave number k with collision 
frequency v = 0.05. The curves are the evolving of the square root of Eh in time using 
logarithm scale. The slope of the line is obtained by the least square fitting. 
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